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We consider a thermofield approach to analyze the evolution of an open quantum system coupled 
to an environment at finite temperature. In this approach, the finite temperature environment 
is exactly mapped onto two virtual environments at zero temperature. These two environments 
are then unitarily transformed into two different chains of oscillators, leading to a one dimensional 
structure that can be numerically studied using tensor network techniques. 


In the past decades, many different techniques have 
been developed to analyze the dynamics of quantum sys¬ 
tems coupled to an environment, i.e. open quantum 
systems (OQS). Some of these are based on deriving a 
master equation, which evolves the reduced density op¬ 
erator of the OQS by tracing out the environment de¬ 
grees of freedom 1HI2], and some others are based on 
the stochastic Schrodinger equations (SSE), evolving the 
OQS’s wave function conditioned by a continuous US] 
or discrete urn stochastic process. Both approaches are 
suitable for weak system-environment couplings, which 
generally lead to a large separation between system and 
environment time scales. Although such a large separa¬ 
tion often occurs in quantum optics, it is not necessarily 
so in other scenarios, such as soft or condensed matter 
systems, or in quantum biology. In these situations, other 
approaches are more appropriate, such as the path inte¬ 
gral Montecarlo 7j, which in some parameter regimes is 
nevertheless hindered by the sign problem, potentially af¬ 
fecting the convergence of the method at relatively short 
times (see for instance ED- 

An alternative is to solve the total system dynamics 
with exact diagonalization methods, which is difficult due 
to the large number of degrees of freedom in the environ¬ 
ment. Hence, a wise selection of the relevant states of 
the full system is of primary importance, and this can be 
done for instance by discarding states with low probabil¬ 
ity, as in the density matrix approach |J] (closely related 
to density matrix renormalization group), or by consid¬ 
ering as relevant only those states generated during the 
evolution, as done in the variational approach mm- 

Similarly, it is possible to perform a unitary trans¬ 
formation of the environment that maps it onto a one 
dimensional structure. The numerical renormalization 
group (NRG) approach [I2HT7] . for instance, is based on 
a (logarithmic) coarse-graining of the continuous envi¬ 
ronment spectral function in energy space. The result¬ 
ing discretized environment can then be mapped onto a 
semi-infinite tight-binding chain [181 with exponentially 
decreasing couplings. As proposed in [MU], the map¬ 
ping can also be performed analytically without a pre¬ 
vious discretization of the environment. Even when the 
couplings do not decay exponentially, it is typically pos¬ 


sible to describe the system dynamics until its decay or 
relaxation time using a truncated chain of finite length. 
The total system can now be modelled as a matrix prod¬ 
uct state (MPS), and it is then possible to use tensor 
network techniques to simulate the unitary evolution of 
the total system [22H2H] . The approach can also deal 
with an environment at finite temperature, using matrix 
product operators mm- 
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FIG. 1. Fig. a) represents the initial problem described 
with 0 of an OQS coupled to a harmonic oscillator reservoir 
at finite temperature. Fig. b) represents the thermofield- 
transformed problem in which the finite temperature of 
the reservoir is encoded in two different reservoirs at zero 
temperature. Fig. c) is the chain representation of the latter. 

In this letter we present a complementary formulation 
of this idea based on the thermofield approach proposed 
in [[2M3U] (see [3T] for a review). In this approach, the 
environmental Hilbert space is mirrored or doubled, and 
then a thermal Bogoliubov transformation is performed. 
As a result, the real environment in a thermal state is 
transformed into two virtual environments in a vacuum 
state (see Fig. 0b), known in the literature as the ther¬ 
mofield vacuum. The expectation value of any operator 
of the real environment in the thermofield vacuum co¬ 
incides with its expectation value in the thermal state. 
The only excitations appearing in the environment will 
be those that are dynamically created through the in¬ 
teraction, and the dynamics of the resulting transformed 
system can be simulated using MPS. 

The thermofield approach has been considered in the 
framework of SSE of OQS (see for instance EJ) , but most 
of its applications are in the context of quantum field 
theory and general relativity [321 S3] ■ 
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Thermofield dynamics- Let us consider an environ¬ 
ment of harmonic oscillators, with annihilation (creation) 
operators b k (b k ) and frequencies ui k , to which the OQS 
couples with strengths g k . The complete Hamiltonian 
can be written as 

Htot = Hg + y ' u! k h k h k T E g k (L% + blL), (1) 

k k 

where Hg is the Hamiltonian of the OQS, and L is the 
coupling operator acting on the OQS Hilbert space. We 
can introduce an auxiliary, decoupled environment, char¬ 
acterized by annihilation (creation) operators c k (c][) and 
write the total Hamiltonian as 

tftot — Htot y \ OJ k C k C k . (2) 

k 

Assuming now that both environments are initially in 
a thermal state at inverse temperature /3, we apply a 
thermal Bogoliubov transformation, 

aik = e~ xG b k e lG = cosh(0 fc )6 fe - sinh(d fe )<4, 
a 2 fc = e~ lG c k e tG = cosh (9 k )c k - sinh(0 fe )&|.. (3) 

Here, G = i Yf, k @k(b k c k ~ c k b k ), with 9 k a function of the 
temperature such that cosh(dfc) = \J 1 + n k , and n k = 
l/(e^“ fc — 1) is the number of excitations in mode k. In 
terms of these new modes, 

fltot = H s J r'^2,uJ k (a\ k a lk - a\ k a 2k ) 
k 

+ gi k {L^ a\ k + a\ k L) + g 2k {La 2 k + a\ k lf), (4) 
k k 

where, g lk = g k cosh(9 k ) and g 2k = g k sinh(0 fe ). The 
thermal vacuum can be written in terms of the vacuum 
for b k , c k modes, |S2o)> as 

\n) = e~ iG \n 0 ). ( 5 ) 

The thermal vacuum can be written in alternative ways 
that further enlighten its physical meaning. Firstly, it 
can be written as |S2) = e~ s / 2 e^ h |STo), with S = 
- J2 k ( b k t>k logsinh 2 (0 fc ) - b k b\ log cosh 2 (0 fc )), which can 
be interpreted as the entropy operator for the physical 
(original) environment |31] , since the thermofield vacuum 
is the state that minimizes the thermodynamic poten¬ 
tial (f2|(-iS + H) |f2). Secondly, up to normalization, 

|0) oc e~P HB / 2 \ J), where 1 1) = I n)b\n) c is a max¬ 

imally entangled state between the real and the auxil¬ 
iary environments, defined in terms of their energy eigen¬ 
states, | n)b, | n) c . The thermal state of the original envi¬ 
ronment is thus pb = Tr aux [|fi)(n|] and it can be approx¬ 
imated by a MPO by evolving the maximally entangled 
state in imaginary time j26 ] l27|. In contrast, the present 
approach is based on directly calculating the dynamics 
of the whole system under the Hamiltonian Q, using 
the thermofield vacuum as (pure) initial state for both 


reservoirs. Although this state is annihilated by a\ k and 
(i 2 fc, the number of physical particles has non-vanishing 
expectation value n k = (fl|&j)&/c|r2) = sinh 2 (0fc). Hence, 
solving the dynamics of the initial problem .0 with an 
initial condition pg ot = pg (g) p*g, with pg the initial 
state of the system, is equivalent to solving the dynam¬ 
ics with Q, but considering pg ot = pg g) |H)(fl| (see 
Fig. 0) and (JTJd) respectively). We have described in 
detail the thermofield transformation for bosonic envi¬ 
ronments, but a similar Bogoliubov transformation can 
be proposed for fermionic reservoirs. In that case [31] 
we have a\ k = e~ lG b k e lG = cos (9 k )b k — sin(0fc)cj), and 
a 2k = e~ lG c k e lG = cos [9 k )c k + sm(9 k )b\,. With this 
transformation, the Hamiltonian 0 is transformed into 

0; . . 

Chain representation- The Hamiltonian (I4J) represents 

an OQS interacting with two independent environments, 
having operators a\ k and a 2k respectively. The whole 
problem can be mapped into a one dimensional structure 
with the schematic form in Fig. 0 )- In general, the 
environment oscillators in Q form a quasi-continuum, so 
that the Hamiltonian can also be written as H = Hg + 
fg 1 dkg(k)(b(k)L^+Lb(k) t)+ u(k)b{k)'b(k). When the 
environment is in a Gaussian state, w(fc) and g(k) enter 
the description of the OQS only through the spectral 
density, J(w). In this situation, one can always choose 
w(fc) = u>ok (with ujo an arbitrary constant that may 
be taken as one), and g(k) = \JJ(uj(k)). Similarly, the 
continuum representation of Q reads 


iZtot = H s + j dkk(a\ k a lk - a\ k a 2k ) 

Jo 

+ / dk[gi k {L^ a\ k + a\ k L) + g 2k {La2 k + a) 2k L^)\ 
Jo 


Thus, the spectral densities are J\{k) = gl(k) = J^ fe (l + 
n(ui(k))) J(w(fc)), and J 2 {k ) = g 2 (k) = n(w(fc))J(w(fc)). 
Then, using the unitary transformation discussed in m 
[ 20 ] , new bosonic operators B n and C n can be defined for 
each reservoir, such that 

O'lk = y ' Ui n {k)B n , a 2k — y ' U 2 n(k)C n: (6) 

n n 

where U jn (k) = gj(k)ir jn (k)/p nj (j = 1,2). Here, 
7 Tj n [k) are rnonic orthogonal polynomials that obey 
fg dkJj{k)'Kj tn (k)-Kj trn {k) = PnjSnm, with p 2 nj = 
fg dkJj(k)n 2 n (k) [20J [2T]. Hence, the proposed trans¬ 
formation is also orthogonal, f dkU* n Uj m = S nm . The 
transformed Hamiltonian can be written as H to t = Hg + 
HB+H ln t, with H mt = gi{L^B( j +BgL)+g 2 (LCo+CgL 11 ), 
with gj = pjo , and 


Hb — y ' {^cr.\^ n B^ l B n OL 2l nC^ l C n 

n—0,--- ,M 

,n+iBl l+1 B n -^P^Ci +1 C n + h.c.), ( 7 ) 
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where the recurrence relation of the polynomials has been 
used, Trj tn+1 (k) = (k - a^ n )it n (k) - with 

7t j n—i = 0. Coefficients a : j >n and j3j >n can be obtained 
with standard numerical routines 1341 . The resulting 
Hamiltonian describes two tight-binding chains to which 
the system is coupled. The thermofield vacuum is also 
annihilated by the new modes B n and C n , so that the 
dynamics of the whole system can be simulated using 
MPS time-evolution methods from an initial state with 
zero occupancy of each of these modes. 

A similar mapping can be applied in the case of a fi¬ 
nite discrete environment by means of a standard tri- 
diagonalization. 

In the following, we present numerical results to illus¬ 
trate this approach in different examples. 

Example 1: A spin in a bosonic field- Let us consider 
a spin 1/2 system coupled to a bosonic environment with 
spectral density given by the Caldeira and Leggett model 
[551155] , 


J(u) = r]u s e- U/U % ( 8 ) 


with 0 < s < 1 in the sub-ohmic case, and s > 1 in 
the super-ohmic. Roughly speaking, the constant p gives 
the coupling strength between system and environment. 
The exponential factor in ([8]) provides a smooth cut-off 
for the spectral density, modulated by a frequency cut¬ 
off ui c . This general model provides a good approxima¬ 
tion for spectral densities appearing in many different 
problems, like an impurity in a photonic crystal [37, ■'38]. 
quantum impurity models [35], and solid state devices 
at low temperatures such as superconducting qubits m, 
quantum dots m, and nanomechanical oscillators |42| , 
to name just a few examples. As a first check we consider 
a solvable example, with H$ = \hujJs<x Zl and L = a z in 
0. For an initial state \ipo) = a|0) + 6|1), the expecta¬ 
tion value of any system operator, A, can be analytically 
calculated [45] , 

(A(t)j = e~ 2 ^{\a\ 2 e iu3st + |6| 2 e _i “ st }, (9) 


with tfit = fo^ T fo^ s ^ e l a T( T ~ s)], cnT(t) = 
J2k9k coth ( ! ^) cos (uikt) — i sin (wfcf) . To simulate 
the problem numerically using MPS we need to truncate 
the maximum occupation number of the bosonic modes, 
and the length of the chains corresponding to the trans¬ 
formed environment. We compare the numerical solu¬ 
tion to the exact one for A = a x in Fig. [2] and observe 
very good agreement for all considered spectral densities, 
couplings and temperatures, for a relatively small bond 
dimension and length of each chain, M. 

In the following, we consider a problem that is not 
exactly solvable, by choosing L = cr x , and compare the 
solutions of our method with those corresponding to a 
master equation (ME) up to second order in the system- 



t t 


FIG. 2. Evolution of the mean value of (c x (t)) for a = b = 
1 /x/2, for p = 0.1 and u>s = 0. Upper panels correspond 
to the ohmic model (s = 1), for p = 5 (left) and = 1 
(right panel). The two lower panels correspond respectively 
to a sub-ohmic model with s = 1/2 (left panel), and a super- 
ohmic model with s = 3/2 (right panel), both for ft = 1. 
The exact solution is given by the solid black curves. For 
the MPS, we consider a varying M: magenta, green and blue 
curves correspond to M = 2,10,40 respectively in all plots. 
The last plot includes an orange curve with M = 120. Bond 
dimension is D = 20, and maximum occupation numbers in 
the harmonic oscillator basis is n = 3. 


environment coupling parameter g 


= ~i[Hs,Ps(t)]+ J dTa* 2 (t-T)[L\p s {t)L{T-t)\ 

+ [ dra 2 {t - r)[L t (r - t)p s (t),L] 

Jo 

+ [ drai(t-r)[L(r-t)p s (t),L t ] 

Jo 

+ f dra\{t — t)[L, p s (t)L(r — t) f ] + 0(g 3 ), ( 10 ) 

Jo 


with an(t - t) = J2k9k( n k + l)e iuik{t T \ a 2 (t - r) = 
Ex9k n ke iUk{t ~ T) , and L(t) = e iHst Le~ iHst . To derive 
this equation, the Born approximation has also been as¬ 
sumed. This ME neglects the system-environment cor¬ 
relations, and considers that the latter remains in the 
thermal equilibrium state pb during the interaction, so 
that ptot{t) ~ p s (t) <g> Pb- 

As shown in Fig. i> , the ME and the MPS coincide 
quite reasonably at weak couplings. However, as shown 
in Fig. Q, for stronger couplings the ME does not give 
an accurate description of the dynamics. Indeed, the 
MPS results describe comparatively a much slower decay 
for the two temperature values here considered. Also, the 
computational cost of the MPS in the strong coupling 
regime is much higher than at weak coupling. Neverthe¬ 
less, the difference of the present scheme is that the exci- 
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FIG. 3. Comparison of ME (solid curves) with MPS (trian¬ 
gles) results for (a x (t)) considering y = 0.01, uis = 0.1, and 
M = 100 for both chains. We consider 0 = 10 (blue), 0 = 50 
(orange) and 0 = 1000 (red). The MPS results converge with 
maximum population per oscillator n = 4. 

tations involved in the numerical resolution are just those 
that are dynamically generated due to the interaction with 
the OQS. This is in clear contrast with traditional meth¬ 
ods in which the initial state is thermal, and therefore 
already has a finite initial occupation in the environment 
basis. 



FIG. 4. Evolution of ( a x (t )) for 0 = 10 (upper panel, with 
maximum bond dimension D = 40) and 0 = 50 (lower 
panel, with D = 20) for y = 0.1, cos = 0.1, and M = 100 
for both chains. The solid black curve corresponds to the 
solution for the ME. Considering m the dimension of the 
first two oscillators in the chain, and ri 2 the dimension of 
the following ones, the curves with green squares correspond 
to (m = 5,ri2 = 4), and the ones with blue diamonds to 
(m = 6, ii 2 = 5) (lower panel) and (ni = 7, n 2 = 6) (upper 
panel). The curve with orange triangles in the upper panel 
corresponds to (ni = 8, n 2 = 7). 

Example 2: A quantum dot coupled to an electronic 
reservoir- As noted above, our proposal is valid also 
for fermionic environments. To illustrate this, we con¬ 
sider a quantum dot (QD) coupled to an electronic reser¬ 
voir at a finite temperature with a Hamiltonian H = 
H s + H B + H hy . Here, H s = J2a( Vn °- + f n ^ n s) 
is the Hamiltonian of the quantum dot, which is rep¬ 
resented using the Anderson impurity model with an on¬ 
site Coulomb repulsion U and an on-site energy V. Here, 
the operator n a = d]yd a measures the number of elec¬ 


trons with spin a =t,4- a t the dot. We consider that the 
QD is connected to the reservoir through a hybridization 
term y = —t^2 k . a gk{d\,b k + h.c.), that is a sum of 
bilinear terms wherein d\ a (di a ) creates (annihilates) an 
electron at the dot with spin a and b\{bk) creates (an¬ 
nihilates) an electron with arbitrary spin and momen¬ 
tum k in the reservoir. Hence, the interaction Hamilto¬ 
nian has a similar form as the one in |l]), but redefining 
L = —t^pdc. For simplicity, we have considered that 
both spins a couple equally to the reservoir. The Hamil¬ 
tonian of the environment is H B = ')Zk u, kb\.bk. After 
the thermofield transformation, the former Hamiltonian 
is written in terms of H B = X^ w fc( a ifc a ifc - a\ k a 2 k ), 
and an interaction Hamiltonian of the form Q with cou¬ 
plings gik = -tg k sj\ + A-, and g 2k = -tg k y/Tk, with 
fk = (1 + exp(/?Wfc)) -1 . We consider a spectral density 
of sub-ohmic type, with s = 0.5 in Eq. (|8|. 

Comparing the MPS results to those of ME, as shown 
in Fig. [5j we find initial agreement as expected, but then 
the results start to differ considerably even at relatively 
weak couplings. Due to the limited size of the fermionic 
basis, the MPS converges to the exact result with rela¬ 
tively small computational resources. 



FIG. 5. Evolution of (n-|-) (green) and (nj.) (purple) for the 
ME (solid lines) and the t-DMRG (symbols). We have consid¬ 
ered 0=1 (upper panel), 0 = 10 (lower panel) with U = 0.2, 
V = —17/2, w c = 15, t = 0.01, and a M = 100 oscillators in 
the chain. 

Conclusions and outlook - Based on a thermofield ap¬ 
proach, our formalism allows us to efficiently integrate 
the dynamics of an OQS coupled to a thermal reser¬ 
voir, either bosonic or fermionic, in a pure state formal¬ 
ism, without previously preparing the thermal state with 
imaginary time evolution. The approach is based on per¬ 
forming an analytical (thermal Bogoliubov) transforma¬ 
tion over the (physical) environment and an auxiliary 
one. Provided the thermal state of the original envi¬ 
ronment is known, more concretely, that the quantities 
n k can be analytically or numerically computed, our ap¬ 
proach can be used to solve thermalization problems of 
OQS using only zero-tenrperature (pure state) MPS. 
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